'''
Author: Zhuangyu
Date: 2022-09-28 13:20:27
LastEditTime: 2023-01-03 11:44:28
'''
from casadi import*
import numpy as np
from Parameters import *
from RichardsModel_casadi import *

#from Model import *

import random
p=Loam()

## Creating a function that pre-computes the jacobian

Nr,Nt,Nz,dr,dt,dz,Np,NN,Nx,Nw,Ny,Nv,Hz=circular_parameters()
DeltaT,Nsim=time_parameters()

x=SX.sym('x',Nz)
w=SX.sym('w',Nz)
u=SX.sym('u',1)
Kc=SX.sym('Kc',1)
ET0=SX.sym('Et0',1)
DeltaT=SX.sym('DeltaT',1)


def pre_jacobian(x,u,w,Kc,ET0,DeltaT):
    
    F = ODE_discrete(x,u,w,Kc,ET0,DeltaT)
    F_jacx=jacobian(F,x)
    F_jacu=jacobian(F,u)
    
    H = getOutputs1D(x)
    
    ## Computing the jacobian based on the measurement equation
    h_jacx=jacobian(H,x)
    
    return F_jacx, F_jacu, h_jacx

F_symbol1=pre_jacobian(x,u,w,Kc,ET0,DeltaT)[0]
F_symbol2=pre_jacobian(x,u,w,Kc,ET0,DeltaT)[1]
H_symbol=pre_jacobian(x,u,w,Kc,ET0,DeltaT)[2]

F1=Function('F1',[x,u,w,Kc,ET0,DeltaT],[F_symbol1])
F2=Function('F2',[x,u,w,Kc,ET0,DeltaT],[F_symbol2])
Fy = Function('Fy', [x], [H_symbol])
